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ABSTRACT 

The new approach outhned in Paper I (Spurzem & Giersz 1996) to foUow the individual 
formation and evolution of binaries in an evolving, equal point-mass star cluster is 
extended for the self-consistent treatment of relaxation and close three- and four-body 
encounters for many binaries (typically a few percent of the initial number of stars in 
the cluster) . The distribution of single stars is treated as a conducting gas sphere with 
a standard anisotropic gaseous model. A Monte Carlo technique is used to model the 
motion of binaries, their formation and subsequent hardening by close encounters, and 
their relaxation (dynamical friction) with single stars and other binaries. The results 
are a further approach towards a realistic model of globular clusters with primordial 
binaries without using special hardware. We present, as our main result, the self- 
consistent evolution of a cluster consisting of 300.000 equal point-mass stars, plus 
30.000 equal mass binaries over several hundred half-mass relaxation times, well into 
the phase where most of the binaries have been dissolved and evacuated from the core. 
In a self-consistent model it is the first time that such a realistically large number of 
binaries is evolving in a cluster with an even ten times larger number of single stars. 
Due to the Monte Carlo treatment of the binaries we can at every moment analyze 
their external and internal parameters in the cluster as in an TV-body simulation. 
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1 INTRODUCTION 

Dynamical modelling of globular clusters and other coUi- 
sional stellar systems (like galactic nuclei, rich open clus- 
ters, rich galaxy clusters) still suffers from severe drawbacks. 
They are due partly to the poor understanding of the va- 
lidity of assumptions used in statistical modelling based on 
the Fokker-Planck and other approximations on one hand, 
and due to statistical noise and the impossibility to directly 
model realistic particle numbers with the presently available 
hardware, on the other hand. 

Only recently a detailed comparison of the different 
methods for comparable parameter choices has been tack- 
led (Giersz & Heggie 1994a,b, henceforth GHI, GHII, Giersz 
& Spurzem 1994, henceforth GS, Spurzem & Aarseth 1996, 
Giersz 1996, 1998, Spurzem 1996, Heggie et al. 1998). They 
include theoretical models as the direct numerical solution of 
the orbit-averaged Fokker-Planck equation for isotropic sys- 
tems (Cohn 1980), its 2D generalization to anisotropic mod- 
els (Takahashi 1995, 1996, 1997, Takahashi, Lee & Inagaki 



1997) and rotating axisymmetric clusters (Einsel & Spurzem 
1999), isotropic (Heggie 1984) and anisotropic gaseous mod- 
els (Louis & Spurzem 1991, Spurzem 1994) and direct A''- 
body simulations using standard A'^-body codes (NBODY5: 
Aarseth 1985, Spurzem & Aarseth 1996; NB0DY4: Makino 
& Aarseth 1992, Makino 1996, Aarseth & Heggie 1998, 
NB0DY6 and NBODY6-h-h: Aarseth 1996, 1999, Spurzem 
1999) or Monte Carlo schemes (Giersz 1996, 1998). All the 
cited work, however, only dealt with idealized single-mass 
models. There are very few attempts yet to extend the quan- 
titative comparisons to more realistic star clusters contain- 
ing different mass bins or even a continuous mass spectrum 
(Spurzem & Takahashi 1995, Giersz & Heggie 1997). 

The results could be summarized by saying that in gen- 
eral the Fokker-Planck approximation (small angle two-body 
scattering dominates the global evolution of the system), the 
approximation of heat conduction (its energy transport can 
be treated as heat conduction in a coUisional gas), and the 
statistical binary treatment (model of energy generation by 
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formation and subsequent hardening of thrcc-body binaries 
using simple semi-analytical estimates) all appear to be a 
fairly good description of what happens in AT-body Simula^ 
tions. But there are still two basic drawbacks: (i) all com- 
parisons are so far limited to rather small particle numbers 
(A'^ < 64000) as compared to real particle numbers of globu- 
lar clusters of the order of a few 10® or even up to 10® stars. 
Low-iV models cannot be easily extrapolated to higher N, 
since after core collapse a variety of different processes (close 
encounters, tidal two-body encounters, effects of the finite 
size of the stars) all vary with time scales, which depend on 
different powers of the particle number (see e.g. the scaling 
problem tackled by Aarscth & Hcggic 1998); (ii) during core 
bounce and binary driven post-collapse evolution an individ- 
ual iV-body simulation exhibits stochastic fluctuations, due 
to the stochastic occurrence of superelastic scatterings of 
very hard binaries with field stars and other binaries (three- 
body and four-body encounters, henceforth briefly "3b" and 
"4b" encounters). Although the averaged evolution of the 
system, understood either as a time average (looking for 
long post-collapse times) or as an ensemble average (averag- 
ing statistically independent single AT-body models), is re- 
produced well by the theoretical models based on the above 
assumptions, the individual evolution of a stellar system, 
even with a relatively large particle number, might not be 
exactly matched at any instant. The most recent coUaborar 
tive experiment in this area (Heggie et al. 1998) gives a good 
overview: all methods do agree fairly well, but variations of 
quantitative results of some 10 or 20 % and some scaling 
problems, which are not exhaustively examined, have to be 
tolerated. 

Such considerations led to the construction of special- 
purpose computers for direct A''-body simulations (Sugi- 
moto et al. 1990, Makino et al. 1997, Makino & Taiji 1998) 
and considerable efforts to improve the highly accurate N- 
body simulation software used (see e.g. Aarseth 1996, 1999, 
Spurzem 1999, but also the alternative approach by KIRA 
mentioned e.g. in Makino & Taiji 1998, McMillan & Hut 
1996, Portegies Zwart et al. 1998). But there is an elegant 
alternative way to generate models of star clusters, which 
correctly reproduce the stochastic features of real star clus- 
ters, but without really integrating all orbits directly as in 
an AT-body simulation. These so-called Monte Carlo mod- 
els were first presented by Henon (1971, 1975, Spitzer 1975) 
and later improved by Stodolkiewicz (1982, 1985, 1986) and 
in further work by Giersz (1996, 1998). The basic idea is to 
have pseudo-particles, whose orbital parameters are given 
in a smooth, self-consistent potential. However, their orbital 
motion is not explicitly followed; to model interactions with 
other particles like two-body relaxation by distant encoun- 
ters or strong interactions between binaries and field stars, 
a position of the particle in its orbit and further free param- 
eters of the individual encounter are picked from an appro- 
priate distribution by using random numbers. For a more 
detailed description see the cited papers and Sect. 2 below. 

In the past it proved to be very delicate to properly tai- 
lor flexible and versatile pure Monte Carlo models. Stochas- 
tic fluctuations of the parameters of A'^-body systems with 
large particle rmmbcrs (say A'^ > 10000) are mainly due 
to binary effects (formation by 3b encounters, superelas- 
tic scatterings with field stars, or, in the case of the pres- 
ence of many primordial binaries, binary-binary encounters. 



see e.g. Heggie & Aarseth 1992, henceforth HA92), not so 
much due to the inherent fluctuations resulting from the dis- 
crete nature of the particle system (which should decay ap- 
proximately with \/N). Therefore the idea appeared only to 
model the binary population by a Monte Carlo technique, 
above a background of single stars, which are treated by 
a standard theoretical model. Talcahashi & Inagaki (1991) 
published a similar approach for the case of an isotropic 
Fokker-Planck model, but without following individual bi- 
nary's orbits and their relaxation interaction with them- 
selves and other single stars. In an earlier approach Inagaki 
& Hut (1988) simulated the stochastic evolution of a binary 
population including a distribution of binary orbits. How- 
ever, their background single star cluster was only approx- 
imately modelled, and their binaries were assumed to have 
purely radial orbits, and dynamical friction was treated ap- 
proximately (as opposed to a full Monte Carlo model, which 
models explicitly the cumulative effect of many small an- 
gle encounters in order to generate proper relaxation and 
dynamical friction effects). 

Here we want to extend the stochastic treatment of bi- 
naries in an anisotropic geiseous model for the single stars, as 
introduced in Spurzem & Giersz (1996, henceforth Paper I), 
for the self-consistent inclusion of many primordial binaries 
in a globular cluster model (which is yet a simplified model, 
with regard to its neglect of any finite size effect of the stars 
and its assumption that all stars have equal masses). It is 
very clear now from the observational (Hut et al. 1992) as 
well £is the theoretical viewpoint of star formation that glob- 
ular clusters do contain of the order of 5-15 % of binaries, 
which means that at the time of their formation the binary 
fraction had to be even larger, because many soft binaries 
are subject to disruption by close 3b and 4b encounters. Nev- 
ertheless theoretical modelling of the dynamical evolution of 
even idealized large star clusters containing many binaries 
has not advanced very much in the last years. 

Pioneering studies by McMillan, Hut & Makino (1990, 
1991), McMillan & Hut (1994) using a direct N-body in- 
tegrator with regularization techniques for close binaries 
(NBODY5: Aarseth 1985) used some 1000 particles of equal 
mass only (binary membership of the order of 10-20 %). A 
similar work by HA92 provided in some more detail went up 
to N = 2000, with less than 10 % of binaries. The results of 
such models was that the details of the initial distribution 
of primordial binaries, at least for such small total particle 
numbers, strongly influences the global dynamics of the star 
cluster, which is quite understandable, since usually the bi- 
naries' total binding energy is much larger than the binding 
energy of the entire star cluster. 

On the liigh-A^ branch Fokker-Planck modellers tried to 
incorporate appropriate average cross-section for 3b and 4b 
encounters in their simulations, as was done so successfully 
for the 3b case (Lee et al. 1991, Giersz & Spurzem 1994). In 
a pioneering study Gao et al. (1991, henceforth GGCM91) 
published the first self-consistent model of an N=330.000 
star cluster with as many as 30.000 binaries and showed how 
gravothermal collapse (and subsequent oscillations) were de- 
layed by a large factor in time due to previous "burning" of 
primordial binaries. In their model the close binary encoun- 
ters, however, had to be treated in a very approximate way, 
as we will discuss later. 

Unfortunately the direct AT-body modelling of a case 
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like the one in the GGCM91 paper is hardly possible today 
even with the fastest special purpose computers, because 
both the scaling of the computational time for the general N- 
body problem is prohibitive and as well as the regularization 
of the close binaries downgrades the performance very much. 
In this paper we use the hybrid Monte Carlo code presented 
in Paper I, which combines an anisotropic gaseous model 
for the single star component with a Monte Carlo stochastic 
treatment of many individual binaries. In our largest models 
wo are able to follow 300.000 stars and .'30.000 binaries as in 
the paper by GGCM91. Our model is entirely self-consistent 
in the way that both the binaries and single stars are sub- 
ject to their own and the other component's gravitational 
field. The binaries arc treated individually with their orbital 
and internal parameters (eccentricity, binding energy, semi- 
major axis). Relaxation of the binaries with themselves and 
the single stars as well as close 3b and 4b encounters are 
treated in a direct, but stochastic way using approximate 
cross sections for the latter and the standard Monte Carlo 
procedure (Henon 1971, Stodolkiewicz 1982, 1986, Giersz 
1998) for the relaxation. We use the same probability dis- 
tribution for the hardening of the harder binary by close 4b 
encounters as GGCM91 (originating from Mikkola 1983a, 
b, 1984a, b); different to them, and more realistically, wc 
compute the probability for a close 4b encounter to occur in 
a self-consistent way from the local density distribution of 
our binaries. We can also in detail follow the evolution and 
compile balances of all the reaction products in any close 
encounter (singles and binaries). 

At this place some more sophisticated models of the 
evolution of a large binary population should not go unmen- 
tioned, i.e. the Hut et al. (1992) model with a very detailed 
study of finite-size effects and using unequal masses, simi- 
larly Sigurdsson & Phinney (1995), but both using a static 
background of single stars. Another interesting subject re- 
lated to this study is the question of the evolution of binary 
parameters in young star forming clusters in the galaxy, 
which was studied by Kroupa (1995) by using NB0DY5 
for direct A'^-body simulations. Our model will be soon fur- 
ther developed to incorporate multi-mass systems, which 
will render it applicable for all such purposes. 

In the next section we will describe only those features 
of the model which have been added as compared to Paper 
I. For the treatment of the ordinary relaxation process and 
the close 3b encounters we refer to Sect. 2 of Paper I. 



2 SELF-CONSISTENT MONTE CARLO 
TREATMENT OF BINARIES 

2.1 General remarks 

In Paper I stochastic binaries were followed in their motion 
and interactions (two-body relaxation or dynamical friction, 
and close 3b encounters) with the single stars. The only feed- 
back to the single stars, which were treated by the gaseous 
model, was the kick heating due to superelastic 3b scatter- 
ings. For Fokker-Planck models Takahashi & Inagaki (1991) 
published a similar approach but without any motion of the 
binaries in the system. In an earlier conference proceed- 
ings Inagaki & Hut (1988) reported a stochastic model of 
a binary population including binary orbits. However, their 



background model cluster was only approximately modelled, 
the binaries were assumed to have purely radial orbits, and 
dynamical friction was treated approximately (as opposed to 
a full Monte Carlo model, which models explicitly the cu- 
mulative effect of many small angle encounters in order to 
generate proper relaxation and dynamical friction effects). 
To our knowledge their approach was never continued nor 
published elsewhere. 

In order to improve our models to make them applicable 
for systems with many (primordial) binaries one has to cope 
with a number of complications. First, the energy balance 
is changed due to 4b (binary-binary) close encounters, and 
binary-binary small angle relaxation has to be taken into 
account according to Hcnon's Monte Carlo method (Henon 
1971). Second, the energetic feedback effect of binary-single 
star relaxation cannot be neglected any more for the en- 
ergy balance of the single stars. Third, binaries move in the 
system due to relaxation and close encounters, and so they 
change the total mass and potential of the entire system. 
Another complication occurs due to the losses and gains of 
single stars by binary formation and disruption of binaries in 
close 4b encounters. All such adjustments generate a feed- 
back via the potential gradient in Euler's equation to the 
single stars. The first points were relatively easy to imple- 
ment: we adopted the approximate 4b cross sections given 
by GGCM91 (1991) as a working model, to be tested and 
improved later; the treatment of binary-binary relaxation is 
similar to Giersz (1998), except that the local binary den- 
sity is estimated here from the smoothed out binary mass 
(see below) via Poisson's equation. The second and third 
point, however, appeared to be extremely difficult to cope 
with. The problems were due to the very different nature 
of the gaseous and Monte Carlo treatment for the two com- 
ponents. The gaseous model approximates density and po- 
tential by smooth functions using Newton's method to find 
the iterative solution at the next time step. In contrast to 
this the stochastic binaries may experience rather quick and 
unsteady changes of their positions, velocities, and binding 
energies. Incorporating such changes in a naive way into the 
gaseous model (just adding up energy and mass balances 
per time step and applying them as a local source or sink 
term in the continuity and energy equations) yielded catas- 
trophic failure. Thus we had to invent a number of measures 
to soften the transitions and links between the two systems 
without spoiling the physically correct treatment of the com- 
bined system and still preserve the advantages of the hybrid 
method. The physical justification is that the gaseous model 
is to be interpreted as an ensemble average over the evolu- 
tion of many independent individual star clusters. So, when 
applying mass or energy changes we have to properly aver- 
age out stochastic variations of the binaries. Our resulting 
code is technically fine-tuned and, frankly speaking, rather 
fragile according to our experience. However, with this pa- 
per we want to publish results obtained with a final code 
version, which we believe to produce physically correct re- 
sults, in an efficient way, and without any unexplained or 
unphysical procedures. Due to difficulties described above 
we think it is justified to give in the following subsections 
some rather technical information how the implementation 
is done. 
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2.2 Binary mass distribution, relaxation, and 
heating processes 

Our main idea to cope with the mentioned problems is to 
include the mass of each binary in the gaseous model as 
a smooth orbital mass distribution, computed by using the 
orbit parameters of the binary, which are known from the 
Monte Carlo model. For the stochastic treatment, the posi- 
tion of a binary is still determined by the appropriate picking 
procedure (Giersz 1998), and if we plot our results we give 
such positions of binaries. However, in Euler's equation in 
the gaseous model a smoothed out mass is used. Thus lo- 
cal potential fluctuations are reduced, since they now only 
reflect the change of the orbital parameters of the binary, 
not the stochastically picked individual binary position. So 
from the viewpoint of the gaseous model a binary appears 
as a smooth mass distribution, extending from its minimum 
to maximum orbital radius in the cluster. The mass distri- 
bution in the orbit is weighted according to dr/vr- Thus the 
total binary mass Mb{r) contained in a sphere of radius r, 
which enters into the force equation of the gaseous model is 
obtained from 



Mf,(r) 



E 



(1) 



where j is an arbitrary index counting the binaries and 
Adbj (r) is computed by integrating the following equations: 

Mbj{r) =0 for r < rmin{Ej , Jj) 

Mbj(r) = mtj for r > rmax(-Ej, J^) (2) 
with the total mass of the j-th binary mtj; for rmin < r < 



Mb At) 



rubj 



dr 

Vr 



with the orbital period 



dr 



Note that 
2 2 

Vr = 

rubj 



nibjr^ 



(3) 



(4) 



(5) 



depends itself on r. Since a standard numerical integration of 
the integrals in Eqs. ^ and ^ is very difficult near the turning 
points of the orbit (v^ — + 0), we use Henon's Monte Carlo 
procedure (Henon 1971) to pick between 50 to 500 radial 
positions in the orbit (depending on how eccentric it is), 
compute Vr for each radius, and approximate the integrals in 
Eqs. ^ and ^ by a summation over these stochastically picked 
radii. The gravitational potential of the binary component is 
then computed in a standard way from the mass distribution 
using the approximation of spherical symmetry ($(r) 
for r —> oo). A slight inconsistency enters here, because the 
new potential is not yet known in Eq. |5| for the computation 
of the new binary mass distribution (see Giersz 1998). 

One more step is necessary to make the gaseous model 
compatible with many stochastic binaries. If its time step 
in high density phases becomes as small as some fraction of 



the central relaxation time, close encounters and relaxation 
encounters with sufficiently big deflection angles cause sub- 
stantial heating and cooling in the stellar core. Note that the 
Monte Carlo model simulates the cumulative effect of small 
angle gravitative encounters by representative encounters, 
whose deflection angle does not need to be very small. Ap- 
plying the energetic effect of such encounters immediately 
in the gaseous model component for the single stars would 
disrupt the stability of the system, because the time scale of 
the changes becomes much shorter than the local relaxation 
time. Therefore we define a reservoir of energy, which is ap- 
plied to the single stars with a maximum rate of -Ec/irx,c, 
where Ec, trx,c are the core energy and the central relax- 
ation time, which is smoothly distributed in the inner 20% 
of the total mass of the cluster, with a power-law cutoff. 
We interpret this as a simplified model of how the sub- or 
suprathermal reaction products themselves relax with the 
other cluster members. For the sink and source terms of the 
mass we use a similar procedure. With these measures the 
gaseous model code is able to cope with the induced mass 
and energy variations of an arbitrary number of binaries. 

Heating terms contain contributions from the following 
processes: 

• energy feedback of relaxation encounters between sin- 
gles and binaries to the single stars. It is easily accomplished, 
since for each relaxation event in the standard Monte Carlo 
procedure we know the deflection angle and the energies 
of the reaction products after the encounter. It is checked 
whether the single star can escape as a result of its relax- 
ation encounter with the binary. If this is not the case the 
energy it gained is added as a heating (or cooling) source to 
the gaseous model equations. The relaxation procedure it- 
self is the same as described in Paper I, but now the energy 
feedback to single stars is not considered negligible. 

• heating due to close 3b encounters. This is the classical 
heating by three-body binaries, which was first implemented 
in a steady approximation by Heggie (1984) and Bettwieser 
& Sugimoto (1984). Its Monte Carlo realization depends on 
an assumed probability for an encounter to occur and a prob- 
ability distribution of AEb, the change of binary binding en- 
ergy (see Spitzer 1987, and detailed description in Paper I). 
Again, the heating is only applied if the single star reaction 
product of a close encounter does not escape. The treatment 
of this process is equal to Paper I. 

• heating due to close 4b encounters. With many bina- 
ries in the system, the probability of close binary-binary 
encounters cannot be neglected any more. The problem is 
very difficult to solve. However, as in the case of 3b encoun- 
ters, for hard binaries a typical event may be defined (at 
least in the case where all single stars have equal mass). 
The typical close encounter between two hard binaries gen- 
erates a "hardening" (gain in binding energy) of the harder 
binary of the two and a disruption of the other one into 
single stars. Depending on the initial binding energies the 
single stars created may escape. If not they contribute to 
the heating term described above. To determine the proba- 
bility for a close 4b encounter we sort all binaries radially 
according to their most recently picked position in their or- 
bits and group them into pairs. For each pair the probability 
is computed and compared to a random number, insofar this 
is again the standard Monte Carlo treatment (Giersz 1998). 
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The probability Pbb for a close encounter to occur is deter 
mined, however, here by using an ncrv-ansatz: 

2Gm6 



Pbb = mnRlvTei • ( 1 + „ 



(6) 



where nt = pb/mb is the local binary number density (de- 
termined from the smoothed out binary mass via a spher- 
ically symmetric Poisson's equation), the binary mass, 
Vvci the relative velocity of the two selected binaries, and 
Rb is 2.5 times the semi-major axis of the harder binary 
(Stodolkiewicz 1986). In his paper a different ansatz was 
given for Pbb, but we checked that by order of magnitude it 
provides the same results. By comparing Pbb with a random 
number it is decided whether an interaction takes place. If 
so, we use the probability distribution given by GGCM91 
(taken from earlier work by Mikkola 1983 a, b, 1984 a, b) 



G(,) = f,(l + I,^) 



-11/4 



(7) 



which gives the fractional translational or recoil energy Er 
released in the encounter, defined by 



y ■■ 



Eb + Ei 



(8) 



(Eb, E'b being the two binding energies). The increase in 
binding energy of the harder binary is thus 



AE = y{Eb + E'b) + E'b 



(9) 



since the softer binary has to be disrupted. One quarter of 
the recoil energy Er is assumed to be carried off by the re- 
maining binary, three quarter are distributed randomly be- 
tween the two stars. Such a procedure is very approximate, 
and it neglects certain other channels of binary-binary reac- 
tions, such as the formation of stable hierarchies with three 
or more particles (compare Aarseth 1999) . While our Monte 
Carlo model in future work will be capable to include more 
details of 4b encounters, for the model presented here we 
decide to use the same reaction processes as GGCM91, for 
a comparison of our results with theirs. 

We use the Monte Carlo rejection technique to pick ran- 
dom values from a distribution f{z) = z^^'^; in case a value 
is accepted we determine 



V ■■ 



2 1 
7~ 



(10) 



where z is the original random number. This provides the 
above probability distribution G(y). By checking the trans- 
lational energies of the reaction products (one binary and 
two singles) we determine whether they escape or not. If 
the single staxs do not escape, their translational energy is 
added to the heating source for the single star component, 
and their mass is added to the single stars as well. Again, 
as described above for the energy, the mass is taken away 
from or added to the single star component smoothly in a re- 
gion defined by the 20% Lagragian radius with a power-law 
cutoff. 



2.3 Adjusting the masses and energies 

The following processes change the masses contained in the 
single or binary components and are properly accumulated 



and taken from or given to the single stars in form of a source 
or sink term in the continuity equation: 

• formation of three-body binaries; this process is treated 
as in Paper I, but now the mass loss for the single stars is 
properly accounted for as a sink term in the gaseous model 
equations; 

• generation of two single stars as a reaction product from 
a 4b encounter, which do not escape; this is accounted for 
in a mass source term in the continuity equation; 

• mass loss of singles and binaries by relaxation, close 3b 

and 4b encounters. This is mass loss to the entire system 
which causes its total binding energy to decrease. 

• finally, there is an ongoing motion of binaries in the 
system, due to relaxation with each other and with single 
stars (dynamical friction) and due to 3b and 4b encounter 
kicks. 

To check total energy conservation we use a balance 
equation with respect to the total energy Stot of the single 
stars only: 



2Stot = + m(*<. + 2*6) 



(11) 



where $s, $6 are the local potentials from singles and bi- 
naries, respectively, and is the 3D velocity dispersion of 
the single stars. Note that the factor of 2 in front of the bi- 
nary potential here occurs, because the binaries provide an 
external potential and do not belong to the self-gravitating 
single star system. To simplify the following discussion we 
consider the above as an energy balance for one discrete ra- 
dial shell in our model. To get the total energy this has to 
be integrated (or in the discrete numerical model summed) 
over the volume. In this sense n = pAV is the mass of single 
stars contained in one radial shell. Any changes in our model 
cause a change 

25Etot = v'^Sn + 6n{^, + 2*6) + iiSv" + ii{5^s + 2<5*6).(12) 

Checking the total energy of the single stars in our model 
means book-keeping of all processes which create changes in 
one of the quantities of the above equation. 

• first, there is a change of $5 due to relaxation, close 3b 

and 4b kicks of the binaries, binary escape, and escape of 
single stars which is caused by a binary-binary encounter. 
All such processes lead to changes in the orbital mass distri- 
bution of the binaries, but not to direct changes of the state 
of the single stars; thus we have 



25Eu 



2fiS^b- 



(13) 



If however, single stars escape from the system which be- 
longed as single stars to the cluster before (this is only pos- 
sible by close 3b encounters between a binary and a single 
star) we have an energy change of 



2<5£'tot = 2/i5*, 



(14) 



because the adjustment of the self-gravitating system pro- 
ceeds such that $s5m = nS^s- 

• second, there is mass loss of single stars due to recoil 
ejections in close 3b encounters. From our numerical exper- 
iments we find that in the case of such an adjustment the 
system reacts with Sv^ <C 1 and ^sSp, = —nS^s, i e. the 
thermal structure is not changed, and the change in the 
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mass distribution is fully represented by the change of po- 
tential energy. Neg lecting the first term and using the second 
identity in Eq. yields: 

25Et,t = + Sfi{^s + 2$6) + nS^b. (15) 

• finally, there are two processes, which shift mass be- 
tween singles and binaries, but do not change the total mass 
distribution, i.e. $s + 2$;, remains constant. They are for- 
mation of binaries by 3b encounters and disruption of bi- 
naries by 4b encounters leading to two non-escaping single 
stars. Here again the above equation seems to be the proper 
ansatz. However, we observe that in such a case the adjust- 
ment of the gaseous model proceeds in a way that 

26Etot = S^iv^ + Sfi{^s + 2*6) + ^fiS^t. (16) 

This follows from the fact that the mass given to or 
taken from the single stars is not localized but distributed 
smoothly inside the central 20% of total mass, which intro- 
duces an additional shift in both the binaries' and singles' 
potential. Again from numerical experiments we find a good 
approximation of the identity 5$s/2 = —S^b/3. 

The reader should keep in mind, however, that all the 
described rather complicated details about the energy bal- 
ance for the different processes do only matter for the pur- 
pose of checking the conservation of energy in the entire cal- 
culation. They have no influence on the physical events and 
progress of the model itself. But they were essential for us in 
the process of debugging the method and finding proper de- 
scriptions for all processes. In one of our models (Gao's run) 
the total binding energy of the binaries exceeds the binding 
energy of the system by several orders of magnitude, and 
this again is some orders of magnitude larger than the total 
energy error we find after a very long integration time. This 
is only possible if our procedure is correct. 

2.4 Final remarks 

Having described the zoo of different processes and balances 
in the previous subsections should illustrate to the reader, 
that in a stochastic binary model like ours we have a detailed 
and individual book-keeping of all these processes. It ensures 
that we have practically the same amount of information as 
in a large A^'-body simulation with primordial binaries. The 
underlying additional approximations here are just spherical 
symmetry, no relaxation processes other than standard two- 
body, and the cross sections for close 3b and 4b encounters. 
The latter will be improved to a detailed numerical three- 
and four-body integration in the future. A''-body simulations 
with primordial binaries have been published by McMillan, 
Hut & Makino (1990) and HA92. In no case was the binary 
number larger than 150, and the total star number did not 
exceed 2500. We use the models of HA92 as a template for 
our models, to check whether there are differences and how 
significant they are. They (around 1992) needed 2000 h CPU 
time on a 10 Mflop workstation, hence the problem requires 
72 Tflop. So it is a major computational job still today, 
though not a question of months. But our real interest stems 
from globular clusters, which may start with initial data like 
the ones used in the Fokker-Planck models of GGCM91, 
using 300.000 single stars and 30.000 primordial binaries. 



It is questionable, whether this job can be done in a di- 
rect simulation even by using a Petafiops computer. We see 
the justification in our efforts to create stochastic and Monte 
Carlo models of star clusters with many binaries in the fact, 
that our model is with a computational effort of two weeks 
on a Pentium II PC, able to provide a full model for the 
GGCM91 case. In excess it yields a wealth of more details 
of all binary related processes and the full stochasticity of 
the binaries as it could only be deduced otherwise from an 
enormous N-body simulation. Note that some of the approx- 
imation used in our model will be released in the near future 
without fundamental difficulties, such as the cross sections 
for 3b and 4b interactions. On the contrary there are as- 
sumptions, which are rather fundamental for our model, and 
not straightforward to remove. They are spherical symmetry 
and the assumption that only two-body relaxation by small 
angle encounters is important. Such a treatment of relax- 
ation may be justified, because collective processes are less 
important in large A'^-body systems. Regarding the assump- 
tion of spherical symmetry, all Monte Carlo models so far 
employed it, because it makes the potential computations 
much quicker and it always leads to purely deterministic 
orbits. This assumption is a much more severe restriction 
in the moment, because a treatment of axisymmetric sys- 
tems would require a not straightforward reformulation of 
all Monte Carlo methods (similar as in the case of Fokker- 
Planck models, Einsel & Spurzem 1999). 

One final note: the reader making comparisons between 
the results presented in the following sections as well as 
resuhs by HA92, Hut, McMillan & Romani (1992), and 
GGCM91 should keep in mind that all results (except for 
the Fokker-Planck models of GGCM91) are individual rep- 
resentations of a statistical ensemble of solutions; all energy 
generation events for example occur randomly and the out- 
come of such close 3b and 4b encounters is treated here 
purely stochastic. So one should not expect a perfect match 
between all models, but an agreement of global and aver- 
aged parameters of the simulated systems. By using differ- 
ent random number sequences for our models, however, one 
can generate a statistical spread of the results, which is com- 
parable to the differences between Monte Carlo models here 
and direct A-body results in HA92. 

In the following we present different set's of test models 
by applying our new method, which we call Monte Carlo 
runs and Heggie's runs, due to the aim to compare with 
corresponding other works. Finally, we model a star cluster 
with 300.000 single stars and 30.000 binaries, called Gao's 
runs, which has no direct counterpart in the direct simu- 
lations, so we have to rely on the capabilities of our code 
and the comparison with the more approximate GGCM91 
models. 



3 RESULTS AND MODELS 
3.1 Method 

The dynamical equations of the gaseous model were dis- 
cretized on a logarithmically equidistant mesh of 200 points 
ranging over eight orders of magnitude and solved by an im- 
plicit Henyey-Newton-Raphson scheme (see GS). More de- 
tails on the Monte Carlo description of our binaries can be 



© 0000 RAS, MNRAS 000, 000-000 



A stochastic Monte Carlo approach 7 




1SD0D 1 GODO 1 7000 1 8000 1 9000 20000 21000 22000 23000 




1E000 1G0O0 17000 



iOOO 1B000 20000 21000 22000 23000 



Figure 1. Model MC5, central density, total heating of the sys- 
tem, and heating due to 3b and 4b close encounters only, as a 
function of time. 
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Figure 2. Model MC5Q, as Fig. 



Figure 3. Model MC5, central density plus crosses for each close 
3b and 4b encounter as a function of time. The abszissa of the 
crosses is chosen according to IkT/Ef,, in order to illustrate the 
correspondence between events and fluctuations in the density 
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found in Paper I, the previous section and the cited htera- 
ture. All results presented in the paper are (if not explicitly 
stated otherwise) given in standard units, i.e. G — M — 1 
and E = -1/4 (Heggie & Mathieu 1986), where M and 
E denote the initial total mass and translational energy of 
all members of the cluster (including single stars and bina- 
ries). The total energy used for this normalization does not 
include the internal (binding) energy of the binaries. 

3.2 Monte Carlo Runs 

To check the reliability of the stochastic Monte Carlo model 
we studied models of a system consisting of 10^ single stars 
without any primordial binaries. Binaries in these runs are 
created only due to the dynamical 3b interactions between 
single stars. In one run, labelled MC5, binaries only harden 
due to interactions with single field stars, but in another run, 
labelled MC5Q, they also interact between themselves. Both 
runs are compared to full Monte Carlo models discussed by 
Giersz (1998). 

The evolution of the central density together with the 
different sources of heating of the gaseous system are pre- 
sented in Figs. |l|and| for models MC5 and MC5Q, respec- 
tively. The curve labelled "3b -I- 4b" shows the heat trans- 
ferred to the single stars (gaseous component) due to binary 



Figure 4. Model MC5Q, as Fig. 



hardening (connected with their interactions with field sin- 
gle stars and other binaries). The curve labelled "Total" 
shows the total heating of the gaseous system including "3b 
+ 4b" contributions (see above) plus an additional heating 
or cooling due to the motion of binaries in the system. Note 
that binaries from the point of view of the gaseous model 
form an external system. Their motion caused by small angle 
two-body interactions with single stars (relaxation process) 
induces changes of the total potential felt by the single stars. 
The total heating is always larger than that connected only 
with the strong binary interactions ("3b -I- 4b"), because 
the binaries preferentially lose their transational energy in 
the relaxation process and sink to the centre of the system. 
Therefore the relaxation interaction with the binaries is al- 
ways, in an averaged sense, a heating (the "Total" curve 
is always above the pure "3b -I- 4b" curve). Note that in- 
dividual relaxation encounters between binaries and single 
stars may lead to cooling also. This is in agreement with the 
expectation from the physical picture of mass segregation. 

Large scale quasi-periodic oscillations are the most pro- 
nounced features of these two figures. They span several or- 
ders of magnitude in central density and have strongly vary- 
ing periods. A characteristic feature of these oscillations is 
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the lack of binary energy generation during the long phases 
of maximum expansion. It is widely accepted that this to- 
gether with a temperature inversion (during these phases) 
are the most important signatures of gravothermal oscil- 
lations (Bettwieser & Sugimoto 1984, McMillan & Engle 
1996). 

In Figs. ^ and ^ crosses show the events of strong inter- 
actions between binaries and field stars and other binaries, 
together with the time evolution of the central density. It is 
clear that most of these interactions take place during the 
phases of maximum density, consistent with the picture of 
gravothermal oscillations. Here our runs also agree with re- 
sults concerning the oscillations observed by Takahashi & In- 
agaki (1991) for a stochastic Fokker-Planck model (stochas- 
tic binary formation and energy generation of binaries, but 
no binary orbits and relaxation interactions with themselves 
and with single stars considered), by Makino (1996) for A''- 
body simulations and by Giersz (1998) for full Monte Carlo 
simulations. Note that maxima of the central density ob- 
served in our model are a few orders of magnitude larger 
than in the full Monte Carlo model. This is connected with 
the way in which the central density is computed in the lat- 
ter. It is estimated from the position and masses of a few 
innermost stars (for details see Giersz 1998), which can lead 
to an underestimation of the density. On the other hand, 
the high spatial resolution of the gaseous model in the core 
(where the innermost shell could consist of a mass of only a 
small fraction of a single star) , and its extremely short time 
step during maximum collapse, both lead to very high cen- 
tral densities. It is worth mentioning that the refrigeration 
cycles (characteristic feature of gravothermal oscillations for 
continuum models already shown in Bettwieser & Sugimoto 
1984) are also present in our models. Though they are very 
noisy, at least four distinct cycles can be separated, which 
correspond to four different oscillation periods. Such periods 
occur from period doublings as they were studied when in- 
creasing A'' in detailed Fokker-Planck and gaseous models of 
gravothermal oscillations (Heggie & Ramamani 1989, Bree- 
den et al. 1994, Spurzem 1994). The shape and the size of 
the loops in the refrigeration cycle are consistent with those 
in the pure Monte Caxlo runs, though the latter showed 
even more fluctuations. It is remarkable, that the stochastic 
events occurring in both Monte Carlo models do not disturb 
much the path of the entire system in phase space, as it is 
provided by the continuum models. 

In Figs. I for model MC5 and Figs. || for model MC5Q 
we provide snapshots of the binary distributions in energy- 
radius space of the cluster for four difi'erent time sets (see 
labels in the figures). The crosses shown in the figures rep- 
resent all binaries which have interactions during the stated 
times (the number of crosses depend on the state of the sys- 
tem - large density means more interactions, small density 
less interactions) . For model MC5 a bimodal distribution of 
binaries can be clearly seen - there are binaries with high 
binding energy far from the core and binaries in the vicinity 
of the core with a wide distribution of binding energies. In 
the course of evolution the build up of a reservoir of binaries 
in the outer halo (binaries in so-called parking orbits, as 
noted already by Hut, McMillan & Romani 1992) is clearly 
visible, and will be even more pronounced in models with 
many primordial binaries presented here in subsection 3.4. 
Such binaries have orbits extending very far out into the 




Figure 7. Comparison of two-component gaseous model with 
stochastic Monte Carlo binary model; central densities of single 
stars and binaries as a function of time, fluctuating curves belong 
to stochastic model, see further explanation in the text. 



halo, and do not enter the core (where strong interactions 
preferentially take place). Binaries on parking orbits can be 
seen in the snapshot as vertically aligned sequences of cross 
symbols at the same constant binary binding energy. This 
occurs because in the Monte Carlo procedure the binaries 
in a nearly invariant orbit are repeatedly picked in differ- 
ent positions. Newly formed binaries appear in the bottom 
left corner of the figures and then in the course of evolution 
(interactions with single stars) they are migrating in the di- 
rection to the right and upwards. Finally they are ejected 
in a single strong interaction or put in a parking orbit in 
a less strong interaction. Note the rather large number of 
binaries in parking orbits. In the full Monte Carlo models 
a much smaller number of such binaries is observed (only 
a few). The reason for this discrepancy is not clear, and 
may be connected with the already mentioned observation 
of higher central densities (leading to larger binary forma- 
tion probability) in the maxima for our model as compared 
to pure Monte Carlo and A'^-body models. For model MC5Q 
the build up of the reservoir is also clearly visible, but with 
smaller number of binaries in parking orbits. Because of the 
presence of 4b encounters, which each destroy one binary, 
the number of binaries in the system drops faster, making 
the strip in the bottom of the figures less populated. Never- 
theless, the bimodal distribution of binaries (parking orbits 
and core and its vicinity) is still well visible. The presence 
of strong 4b interactions in the system creates a wider dis- 
tribution of binaries in energy (note that in each 4b interac- 
tion, while one binary is destroyed, the other is considerably 
hardened) as can be seen for example from comparison of 
the upper right plot of Fig. ^ and Fig. ^ for model MC5 and 
MC5Q, 

respectively. 

In the next two subsections we will discuss the Monte 
Carlo stochastic models with many primordial binaries and 
compare their results with data available in the literature 
(HA92, GGCM91). 
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Figure 5. Model MC5, binary distribution in radius over core radius vs. binding energy in kT plane, for four different times as indicated 
in the plots; note the distinct sequences occurring due to the oscillations of the central density. 
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Figure 8. Model S, core and half-mass radius of the single stars 
as a function of time. 



3.3 Heggie's Runs 

First, we have done one more test model in addition to those 
already discussed in Paper I, checking the correct mass seg- 
regation rate of a two-component system consisting of stars 
of masses mi — 1 and m2 = 2. Hence the mass ratio is 




200 400 600 800 1000 1200 



Figure 9. Model S, Lagrangian radii for the single stars and bi- 
naries (radii containing the innermost 5%, 50% of the single stars, 
and 20%, 50% of the bound binaries remaining in the cluster) as 
a function of time. 

two, which mimics the situation of Heggie's model (HA92, 
see following Sect.) for which single stars and binaries are 
all composed of stars of the same mass. In Fig. |^ we show 
the central density evolution as a function of time for a two- 
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Figure 6. Model MC5Q, binary distribution in radius over core radius vs. binding energy in kT plane, for four diflferent indicated times; 
note the distinct sequences occurring due to the oscillations of the central density. 
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Figure 10. Model S, D, E, and SS; number of bound binaries as 
a function of time. 



Figure 11. Model S, energy balances as a function of time in 
Af-body units. The four different contributions axe described in 
the text. 



component gaseous model (see Spurzem & Takahashi 95, 
Spurzem 1992, the treatment is also equivalent to the same 
kind of gaseous model test runs presented in HA92) and for 
stochastic Monte Carlo model (using a gaseous model only 
for the single stax component). To isolate the mass segrega- 
tion effects only, all close 3b and 4b encounters were artifi- 



cially suppressed in our model. As one can see there is an 
excellent agreement between both models over many orders 
of magnitude in the central binary density. The fluctuations 
of the binary density in the stochastic model are not present 
in the continuous gaseous model. There is also a few percent 
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Figure 13. Model S, D, E, and SS; total number of escaping 
stars in singles and binaries as a function of time. 



difference in the final core collapse time. It is of the same 
order as usual differences between collapse times of either 
different physical models (Fokker-Planck, A'^-body, gas), or 
even the scatter of collapse times between stochastic models 
themselves starting with different realizations of the initial 




Figure 14. Model S, D, E, and SS; total external energy of es- 
caping singles and binaries as a function of time. 
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Figure 15. Model S, D, E, and SS; number of escaping binaries 
as a function of time. 
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Figure 16. Model S, central escape speed as a function of time. 

model (using other random number sequence, see a detailed 
discussion of this effect in GHI). 

Now we discuss a number of simulations which may be 
considered as test cases for our stochastic binary treatment 
in the refined version of the stochastic Monte Carlo method 
(as compared to Paper I). We denote these models by S, 
D, E, and SS, with the same naming as given in Table 1 
of HA92. They are runs with 2500 stars each, containing 75 
primordial binaries with binding energies in the range 2 to 20 
kT (logarithmically equally distributed), where the spatial 
distribution of the binaries is the same as that of the single 
stars (Plummer's model). While this describes the standard 
model S, there are model D (with double the number of 
binaries, but the same binding energy distribution), model 
E (double number of binaries, but stretched from 2 to 200 
kT), and model SS (only 75 binaries as the standard model, 
but binding energy distribution extended from 2 to 2000 
kT). 

All our results match rather closely those of HA92, 
though each run takes only of the order of one hour on a 
Pentium II computer. We do not want to elaborate more 
on the physical interpretation, because this can be referred 
to HA92. Instead we present a number of figures in exactly 
the same manner as there, such that the reader can judge 
how well an agreement can be reached by our Monte Carlo 
type model. Clearly at such low particle numbers we cannot 
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Figure 17. Model S, binary snapshots in a diagram plotting ra- 
dius over initial core radius against binding energy in kT (initial), 
for the first 500 time units. There are more points than actual bi- 
naries because data of several time outputs have been plotted 
together. 
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Figure 18. Model S, as Fig. 0, but for time 1000 < t < 2700. 



expect a complete match of the results. However, we have 
performed several runs for each model with varying random 
number seeds and our general conclusion is, that some fea- 
tures are stable while others exhibit stochastic variations. In 
our following arguments we will always try to discuss which 
differences to the A^'-body models we consider as "real" (in 
the sense that they may still point to deficiencies of our 
stochastic Monte Carlo model) and which we consider as a 
result of the statistical variations at small A*'. 

Fig. ^ shows core and half-mass radii of the single stars 
in our system as a function of time (to be compared with 
Fig. 23 of HA92). It should be noted, that for a compar- 
ison with HA92 here and in all other plots related to the 
core radius, one has to take into account that IIA92 use a 
non-standard definition of the core radius. They use a den- 
sity radius rp, similar to the prescription given in Casertano 
& Hut (1985), because it is operationally well defined for 
A'^-body simulations. In contrast to this the standard core 
radius = 9(T^/(47rGpc), which we use here is in their A''- 
body models more difficult to determine (rc is also some- 
times called the King radius, where CTc, Pc denote the cen- 
tral ID velocity dispersion and central density, respectively). 
As is stated in HA92 for example a Plummer model yields 
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Figure 19. Model SS, as Fig. 
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Figure 20. Model SS, as Fig 



rp — 0.77rc; but during the system's evolution this ratio 
changes. Therefore the core radii measured by HA92 (in 
their Fig. 23) in the late stages have a systematic trend 
to smaller values than those from our models. 

Fig. 1^ shows some Lagrangian radii separately for the 
single stars and binaries (see Fig. 8 of HA92). As in HA92 
the segregation of binaries is very rapid up to about 100 
time units. A more stationary phase follows, characterized 
by continuous binary destruction, and energy generation, 
until most of the binaries are destroyed. For a more detailed 
physical discussion see HA92. The somewhat faster growth 
of the half-mass radius in the gaseous model (as compared 
to the A'^-body result) is a typical feature resulting from 
different outer boundary conditions (see GHI). 

Fig. ^ presents the number of binaries remaining 
bound in the system (see Fig. 9 of HA92). In contrast to 
HA92 after time t x 1500 (model E) or t f=i 1000 (model 
S) the steady binary destruction does not slow down, but 
continues in our models to rather small binary numbers. 
In HA92 thereafter the binary number remains larger. The 
variation of such plots with different random number ini- 
tialization of the system has been explored. We found that 
it can explain to some extent the remaining differences, be- 
cause stochastic variations are big in the late phases when 
there are only relatively few binaries left. The model we are 
using for our Fig. for example, underwent core bounce 
relatively earlier than others, which explains the quicker de- 
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Figure 21. Model S, distribution of binding energies as a function 
of time for all bound binaries. The unit of energy is 2kT, i.e. one 
and half times the initial mean kinetic energy of single stars. Each 
line gives the fraction of these binaries with binding energy less 
than the stated value. 
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Figure 22. Distribution of binding energies of all binaries in the 
model E, see Fig. Ell 




Figure 23. Models S, D, and E; fraction of binaries in the core 
as a function of time. 
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Figure 24. Models S, D, E and SS; core radius as a function of 
time. 



struction of binaries. Selecting another physically equivalent 
initial model can yield diflterent collapse times, and very dif- 
ferent features during the post-collapse phase (see GHI) . We 
want to stress, however, that the global energy budget, pre- 
sented for our model S and D in Figs. and IS m good 
agreement with the features seen in the direct A''-body model 
(see Figs. 11 and 19 of HA92). For these plots we measure 
the total internal (binding) energy of the binaries remain- 
ing bound in the system (E^^^) and of those binaries which 
escaped {E°^f.), and the total external (i.e. potential plus 
translational) energy of all objects (singles and binaries) re- 
maining bound in the system (i?ext) and escaping (£^oxt)- 
Note that i^g^t would be the standard total energy in a sys- 
tem without internal degrees of freedom (hence it starts at 
the canonical value of -0.25). As in HA92 we observe the ini- 
tial phase of an increase of total energy due to hardening of 
the binaries in the system (until t ~ 1000 in Fig, ^l|) , which 
we would like to call the hardening phase, while afterwards 
total energy is gained in part by encounters leading to an 
escape of binaries with high binding energy. The total inter- 
nal energy of bound binaries increases in this phase, which 
means that on average the typical binary in the system be- 
comes less bound. We would like to call this the escaper 
phase. The reader is referred to a much more thorough dis- 
cussion of the physical processes and balances determing 
these phases in HA92. Here we just want to stress that our 
model is able to reproduce all features of the A'^-body model. 

Now we get to Fig. where we have plotted the to- 
tal number of escaping stars, where escaping binaries have 
been counted as two stars (see Fig. 12 of HA92). Compared 
to our results HA92 find much larger numbers, and no dif- 
ference between model S and model D. They conclude that 
the presence of binaries does not enhance the rate of escape 
of single stars very much. In our isolated gaseous models 
we have not yet incorporated single star escapers due to re- 
laxation and tidal fields, which is a subject of future work. 
This means that differences in the escape rates between our 
four models S, SS, D, E should be interpreted only in re- 
lation to the close 3b and 4b encounter processes. Hence 
our model reveals differences regarding the origin of single 
star escapers, which cannot easily be seen in the complete 
A^'-body of HA92, because their models are dominated by 
standard relaxation escapers. In contrast to the number of 
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escapers, the external energy of escapers should in all models 
be dominated by the highly energetic 3b and 4b encounters. 
Consistent with this, our results for the external energy of 
all escaping objects (see Fig. ^ and Fig. 13 of HA92) agree 
much better with the HA92 model. The higher external en- 
ergy of escapers for model E and SS can be attributed to 
the larger maximum binding energy of their binary energy 
distribution, which leads to very high recoil energies in 3b 
and 4b encounters. Also the number of escaping binaries, 
in our Fig. |l5| matches well the HA92 results again (their 
Fig. 14). 

The central escape speed as a function of time, which is 
a direct measure of the central potential, is presented in Fig. 
[l^ , as compared to Fig. 15 of HA92. The A^-body models of 
HA92 reach significantly deeper values of the central poten- 
tial. We think that the reason for this, at least partly, is a 
bias originating from the pairwise potentials of the single 
stars. Hence the potential computed in the direct A'^-body 
model is more prone to clumps of particles, which we do not 
have in our gaseous model representation. Also the max- 
imum is at a different time, in our model at the time of 
maximum binary segregation (roughly at t ~ 200) to the 
centre, in their model at the time of maximum collapse of 
the single stars {t ^ 500). 

In Figs. |l^ to ^ we provide snapshots of the binary 
distributions in energy-radius space of the cluster, for two 
different times, early just in the hardening phase, and late 
well into the escape phase, for model S and SS. In the late 
phases a bimodal binary distribution is observed, soft and 
hard binaries deep in the core, and a group of hard and very 
hard binaries being ejected into the very outer regions of 
the cluster. Very similar behaviour can be seen in Fig. 22 of 
HA92, which gives the corresponding A^-body data. 

Figs. ^ and ^ show the distribution of binding energies 
of the binaries as a function of time for models S and E (see 
Figs. 16 and 20 of HA92); both results again agree rather 
well, except for the bins of the most bound binaries (between 
32 and 128), which are more abundant in our model, due 
to the relatively early creation of new three-body binaries 
in the core. The time at which such first creation of new 
binaries occurs can vary strongly as a function of the random 
number initialization of the model. Note that the units of 
binary binding energy are different by an approximate factor 
of 2 between HA92 and our models; corresponding lines can 
be identified from f = 0. 

To complete our survey of results we present the frac- 
tion of binaries in the core (Fig. and the core radius 
(Fig. 1^) from the four models S, D, E, and SS (see in HA92 
Figs. 17 and 18). A quantitative comparison is difficult due 
to very noisy data of A^-body and stochastic Monte Carlo 
models and due to the different and incompatible core defi- 
nitions in both models. However, the initial increase of the 
binary fraction in the core, and subsequent slow decrease, 
are clearly visible in Fig. G3. In Fig. also the initial fast 
decrease of the core radii for all models is clearly seen, but 
the further evolution, particularly for model D is not fully 
compatible with results of HA92 (see their Fig. 18). Note 
that the above mentioned features, which are in good agree- 
ment, are related to the initial fast mass segregation phase 
of the binaries. 




Figure 25. Evolution of the core and half-mass radii (scaled by 
the scale length of the Plummer model) as a function of time in 
units of initial half-mass relaxation time; this time unit is used in 
all following figures. 
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Figure 26. Mass fraction remaining in single stars and binaries 
as a function of time. 



3.4 Gao's Runs 

Gao's runs are named after the pioneering paper of Gao et 
al. (1991, GGCM91), which we use as a reference paper. 
This is the only paper providing models of live and self- 
consistent single and binary stars in a large number (30.000 
primordial binaries) in a star cluster consisting of 360.000 
stars (including the binary members). Their models were 
two-component Fokker-Planck models, treating the binaries 
as a smooth evolving mass distribution (since all stars have 
the same mass, only two-components are sufficient). The 
initial primordial binary distribution has a binding energy 
range from 3 kT to 400 kT, their spatial distribution is a 
Plummer model realization using the same scaling radius 
as for the single stars. In this subsection we describe, how 
we redo their calculation with our stochastic Monte Carlo 
model, what features can be reproduced, and which ones 
are different. We include the same physics as in GGCM91 
as far as possible and reasonable. This means for example 
we use exactly the same 4b encounter cross-sections that 
they did. However, our model provides more details of the 
binary evolution and of the close encounters. For example, 
about the single star reaction products, which they could not 
provide due to the statistical treatment of their binaries. 
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Figure 27. Ratio of central densities in binaries and singles as 
a function of time. Arrows indicate the times for snapshot plots 
of 2D and 3D distributions of binaries bound to the system. The 
selection of times is explained in the text. 
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Figure 28. Central escape speed as a function of time. Arrows 
as explained in the Fig. [27j 
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Figure 29. Evolution of the core radius of single stars and 1% 
Lagrangian radius of binaries as a function of time. 



There is one technical remark to be made, regarding 
the very important choice of a proper relaxation interval for 
the Monte Carlo binaries (compare standard description in 
Paper I). In contrast to Heggie's runs, here phases occur 
where there are very high densities of single stars and bina- 
ries. Binaries may stay in very elongated orbits, so if we pick 
binaries only at some fraction of their orbital time they are 
being "left behind" in an unphysical way in collapse phases 
of the single stars. In such a situation the results did not 
agree with other models or physical expectations. For those 
binaries newly obtained positions were very far outside, and 
practically no evolution of the binary population took place. 
Therefore we employ (after some experiments) the prescrip- 
tion to use the local relaxation time at two times the core 
radius, divided by ten, as the proper interval for relaxation. 
It would not be harmful to use an even smaller interval, but 
then more and more computational time would be spent to 
relax the binaries. A too small relaxation interval, however, 
would get into conflict with the Monte Carlo method's prin- 
ciples, because the picking of the binaries position at any 
point of the orbit becomes unphysical if done in a time in- 
terval much smaller than the orbital time scale. 

Fig. ^ shows the time evolution of the scaled core and 
half- mass radii of the single stars, to be compared with Fig. 1 
of GGCM91. Like in their model the evolution can be di- 
vided into a first, rather fast, mass segregation phase of the 
binaries (as in the case of Heggie's models), a binary burn- 
ing phase, in which gradually most of the binaries are de- 
stroyed or ejected by close 3b and 4b encounters, and finally 
a standard phase of gravothermal oscillations follows, which 
is dominated by the single stars only with the formation of 
a few new binaries due to 3b encounters. In this phase there 
are only a few hundred primordial binaries left in the out- 
skirts of the cluster, as will become clear later. The model 
continues until about 280 initial half-mass relaxation times 
(trho )• To assess what a heroic effort such a model would 
be in the direct A'^-body approach just note, that this cor- 
responds to 0.8 million A-body time units or some 260.000 
initial half-mass crossing times. We stop the model at about 
280 trho , because nearly all binaries have been destroyed or 
ejected from the core, and the simulation goes on as a single 
star gravothermally oscillating model, where only binaries 
are formed in the high-density phases in the core, as in the 
Monte Carlo runs discussed in subsection 3.2. Different to 
those models, we only have a certain cloud of a few hundred 
binaries in parking orbits very far in the halo, but they do 
not influence the core evolution. Gao's runs reached a sim- 
ilar state already after 90 trho ; but note that the fraction 
of time spent in the different phases (such as oscillations 
or binary burning) seems at a first glance not much differ- 
ent from Gao's case. We will come back to this time scale 
problem below. 

First, we would like to refer the reader to the following 
figures, which show the mass fraction remaining in singles 
and binaries bound to the system (Fig. p^ ), the ratio of 
central densities of binaries over single stars (Fig. [27[ ), the 
central escape speed as a measure of central potential (Fig. 

), and the evolution of the core radius of single stars and 
the 1% Lagrangian radius of the binaries (Fig. ^), all as 
a function of time in units of f^ho • Corresponding figures 
in GGCM91 are Figs. 2 (mass fractions), and 3 (ratio of 
central densities). We will now deduce the scenario of the 
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evolution of our system from these figures, with additional 
information from pictures of the detailed binary distribution 
in energy (over initial kT) and radius (over initial core ra- 
dius), provided for key time points in Figs. ^ and The 
time points are marked by arrows with numbers to 5 in 
Figs. ^ and corresponding times in trho are given in the 
caption of Fig. We will denote the binary snapshots at 
these times in short with snapshot number to 5; they show 
the frequency of binaries (as a 3D surface, with contour lev- 
els projected onto the plane below) over a logarithmically 
equally spaced mesh in scaled energy and position (Fig. , 
or a direct projection of these data (one cross per individual 
binary) in Fig. We deduce the following basic scenario: 

First, up to about 10 trho we have a fast binary mass 
segregation phase, in which the total number of binaries does 
not change much, since the time scale for their destruction 
and ejection is still too long, but as the binary distribution 
flattens out into core and halo, the maximum of the distri- 
bution becomes shallower. The effect can be clearly seen in 
the difference between snapshot and 1, and by the pro- 
nounced maximum of central escape speed reached at time 
point 1 in Fig. |2^. Subsequently there follows a phase of 
rapid binary destruction, as can be seen in Fig. EGI up to 
about 30 trho ; remarkably from the snapshot 2 it can be 
seen that preferrably the low energy binaries (up to a few 
tens of kT) are being destroyed, while the high energy bi- 
naries are not afi'ected much. Here we see the 4b encounters 
being efi'ective; with a high probability the low energy bi- 
naries are undergoing a close encounter with a moderate or 
high energy binary, leading to the destruction of the first 
and hardening of the second binary. Depending on the en- 
ergy of the second binary, escape of single stars only or of 
single stars and the hard binary may happen; if in the be- 
ginning there are still enough binaries of moderate binding 
energy the dominant effect is just heating of the surrounding 
core by 3b and 4b encounters. Therefore we observe a quasi- 
stationary binary burning phase lasting up to about 60 Uw 
, during which the ratio of binary and single star density is 
kept nearly constant in the core (Fig. |27| ) , with a gradually 
expanding core radius of the single stars (Fig. ^) . At about 
60 trho the reservoir of binaries to destroy and heat the core 
is substantially depleted, the binaries have been reduced in 
number by 60%. 

The core starts its first attempt to collapse gravother- 
mally, which is visible as a transient recoUapse of the core 
radius, a shoulder in the binary mass fraction [dM/dt of the 
binaries becomes significantly smaller temporarily), and an 
increase of the ratio of binary to single star density, as a con- 
sequence of mass segregation, all three effects clearly visible 
between 50 and 70 trho in Figs. to However, there 
is still a large enough fraction of the initially 30.000 bina- 
ries present to halt core collapse at a higher level of binary 
to single star density, and cause a second quasi-stationary 
phase until about 150 trho • The process repeats itself, but 
now at a small enough binary number that the single stars 
can start to undergo the first pronounced oscillatory peaks 
(small core radius). Gravothermal oscillations follow, which 
are pronounced at the time points 3 (maximum density), 4 
(expanded stage) and 5 (terminal point of our model). It 
is interesting to note, that the gravothermal oscillations are 
not only visible in the core radius and Lagrangian radii of 
the single stars. They are also visible in the binary density. 



and in the ratio of binary to single star density (see Figs. ^ 
and These features, including amplitude variations over 
many orders of magnitude, multi-peak structure with period 
doublings at the maxima, and long inactive expanded phases 
are clear and the same as observed in GGCM91 and other 
standard models of gravothermal oscillations. From Fig. |30| 
time points 3 and 4 we can see, that in the collapsed phase 
(time point 3) some binaries are very deep in the core, while 
in the expanded phase (time point 4) the core is void of bina- 
ries, no activity taking place. Most interestingly at the final 
time point 5, we have reached a situation where all binaries 
are very far outside in parking orbits. We note, however, that 
from the studies of MC runs (MC5 and MC5Q, subsection 
3.2) the reason, why there are so many binaries in parking 
orbits is not yet clear. In contrast to MC runs, however, here 
most binaries in parking orbits are still primordial and have 
not been created in 3b encounters in high density phases. 
Therefore we believe that the existence of many binaries in 
such orbits is real, though our model may overestimate it. 
In any case, such binaries decouple from the evolution of the 
rest of the system, which will behave more like a single star 
cluster with gravothermal oscillations. Binaries left over in 
the outer halo may be regarded as fossil records of the early 
primordial binary generation. They will be slowly depleted 
by 3b and 4b interactions when they enter the core due to 
relaxation processes. Note that in phases 3 to 5 the scale of 
the binary distribution in Fig. ^ has been amplified by a 
large factor; the absolute frequency (and total number) of 
binaries are one to two orders of magnitude smaller than at 
time points to 2, as can be deduced directly from Fig. ^ 
Figs. ^ and ^ show the evolution of selected La- 
grangian radii for binaries and singles, the number of escap- 
ing stars in binaries, singles, and total (counting binaries as 
two stars), and the total energy budget, in the same way as 
discussed for Heggie's models. Due to their stronger central 
concentration the binaries take part in the gravothermal os- 
cillations with more than 50% of their total mass; at the 
end the final complete removal of binaries from the core can 
be seen in Fig. ^ Figs. ^ and |m| show that the initial bi- 
nary destruction phase is accompanied also by a heavy loss 
of binaries with high binding energy, as can be seen from 
the strong drop of -Ef^^. In the following phase the rate of 
escape becomes slower, but the mechanism is the same, en- 
ergy bound in binaries is carried away by escapers. Note the 
signature of the oscillations in the binary escaper energy at 
late times, and the enormous amount of energy exchanged 
via the binaries, as compared to the initial and final total 
external energy of the star cluster, which starts at the stan- 
dard value of -0.25. Remarkably we lost in total at the end of 
our simulation about 90.000 stars escaped altogether, which 
is one quarter of the initial number of stars. This number 
would be even larger if escapers due to relaxation eff'ects of 
single stars were allowed. Since we lose only some 15.000 bi- 
naries by escape (Fig. ^) , it can be concluded that another 
15.000 binaries have been destroyed by close 4b encounters, 
and, more interestingly, at least 30.000 single stars have es- 
caped which do not originate from one of the primordial 
binaries, assuming that the maximum of all destroyed bina- 
ries led to two single escapers, which is an upper limit only, 
and neglecting possible exchange reactions. In other words, 
per binary destroyed or ejected about three single stars es- 
caped on average. Looking at Fig. 2 of GGCM91 they find 
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Figure 30. Snapshot plots of 3D distributions of bound binaries. Plotted is the frequency of binaries in logarithmically cquiSlistant 
bins of binding energy in kT (initial) and position in units of initial core radius, for six different times in units of the initial half-mass 
relaxation time; top left 0, top right 5.21, middle left 46.51, middle right 228.70, bottom left 246.99, bottom right 278.36. The selection 
of times is the same as indicated by the arrows in Fig. BSl 



an increased number of single stars during the first evolu- 
tionary phase, which cannot be seen in our results. We think 
that this difference is an artifact of the GGCM91 runs, due 
to their rather artificial procedure to select binaries from en- 
ergy bins for close encounters, while in our paper the proper 
encounter probabilities are used. Thereby we destroy less 
quickly the soft binaries and find a higher fraction of es- 
caping single stars from the beginning, so in our results the 
number of bound single stars decreases. Another reason why 
GGCM91 find more bound single stars is that they can not 
take into account the external energy of the destroyed binary 
in a 4b encounter, as a further contribution to the possibly 



escaping single stars in excess of the recoil energy obtained 
by hardening the harder of the two binaries (though the 
latter will in most cases be the dominant contribution). 

Hence, we can explain the time scale difi'erences between 
our and GGCM91 models; we keep a larger number of soft or 
intermediate binaries in the system which are able to supply 
sufficient binary heating, by single stars originating from 4b 
encounters and remaining in the core. These heating reac- 
tions support our long initial binary burning phase, which 
takes place in two phases (as discussed above) until 150 
frho • Opposed to this, the soft and intermediate binaries in 
GGCM91 are quickly destroyed, leading to a much stronger 
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Figure 31. The same data as in Fig. |30| but in a 2D projection of the individual data of each binary (being represented by a cross) 
onto the energy-radius plane, all units as in Fig. 



collapse of the single star and the remaining relatively hard 
binaries, and further acceleration of the recoil and escape 
process. Consistently one could see that the gravothermal 
oscillations in GGCM91 already begin at a fraction of bi- 
naries left of about 30 to 40%, while in our case we have 
a binary depletion to 10% until the oscillations start. Since 
they destroy the weak binaries too quickly this difference in 
evolution between the models occurs, although the time evo- 



lution of the total binary number is not so different. All this 
leads to a time, to begin gravothermal oscillations, which is 
about a factor of 3 shorter than in our models. This is a big 
difference, which however, can be fully explained by a much 
faster binary destruction in their model. We attribute this 
to the rather artificial procedure by which GGCM91 decide 
whether a close 4b encounter takes place. 
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Figure 32. Evolution of Lagrangian radii containing 50% and 
5% of the mass of single stars and 50% and 20% of the mass of 
binaries as a function of time. 
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Figure 34. Energy balance as a function of time. The four dif- 
ferent contributions are described in the text, and are the same 
as used for Figs. IVu and h2l for Heggie's runs. 




... 

Figure 33. Number of escaping single stars, binaries and the 
total number of escapers as a function of time. 

4 CONCLUSIONS AND DISCUSSION 

The new approach outlined in Paper I (Spurzem & Giersz 
1996) to follow, in the Monte Carlo manner, the individual 
formation and evolution of binaries in an evolving point- 
mass cluster was successfully extended to the fully self- 
consistent treatment of relaxation and close three- and four- 
body (abbreviated 3b and 4b) encounters for a substantial 
number (typically a few percent of the initial number of 
stars) of binaries and a realistic total number of stars. We 
use a standard anisotropic gaseous model (Louis & Spurzem 
1991, Spurzem 1994), describing the single stars, and the 
Monte Carlo technique (Giersz 1998) to model, in a stochas- 
tic way, the binary subsystem. The aim of this paper was 
to test the reliability of the new hybrid code by comparing 
its results with the data available in the literature for sin- 
gle mass systems (full Monte Carlo model - Giersz 1998), 
and for systems with substantial primordial binary popula- 
tion (A'^-body models - HA92 and Fokker-Planck models - 
GGCM91). Also we want to show, that our model is able 
to handle large N systems with a large binary number by a 
reasonable effort, but keeps full self-consistency and detailed 
information. It will be possible in the near future, after some 
more real physics is included (such as stellar evolution, finite 
stellar radii, tidal fields, and a mass spectrum), to provide 



unique model data for comparisons with the expected wealth 
of observational data for globular clusters. 

The features of our MC5 and MC5Q runs for isolated, 
single mass systems consisting of 10^ stars with only strong 
3b (MC5) and both 3b and 4b (MC5Q) interactions, re- 
spectively, are in a good agreement with the full Monte 
Carlo model (Giersz 1998). The gravothermal oscillations 
are the most pronounced features of these runs. They show 
very large amplitude in central density, long expanded (in- 
active) phases, and several discrete oscillation frequencies 
which are observed in standard gaseous or Fokker-Planck 
models of gravothermal oscillations (Bettwieser & Sugimoto 
1984, Heggie & Ramamani 1989, Cohn, Hut & Wise 1989, 
Breeden et al. 1994) . Also rapid changes in the phase of the 
oscillations occur, which agree with other stochastic Fokker- 
Planck, A'^-body and Monte Carlo models (Takahashi & In- 
agaki 1991, Makino 1996 and Giersz 1998). As can be ex- 
pected, the expansion phases of the system are powered by 
the themperature inversion in the core, not by the energy 
generated by binaries in interaction with singles and other 
binaries. In the large expansion phases there is practically 
no binary activity, as expected. The binary distribution in 
energy-radius space of the cluster is clearly bimodal. Bina- 
ries with high binding energies and in orbits not entering 
the core and extended far into the halo form one group. 
Binaries with a wide range of binding energies and orbits 
entering the core form the second group. From models with 
artificially suppressed 4b encounters we find that they cause 
a wider distribution of binary binding energies (mainly to- 
wards larger binding energies) and a somewhat less pro- 
nounced bimodal binary distribution in a diagram showing 
the binary distribution in radius and binding energy. Our 
models exhibit a rather large number of binaries in parking 
orbits (orbits which do not reach the core, and have apoc- 
entre very far out in the halo), as compared to full Monte 
Carlo (Giersz 1998) and direct A'-body models (Spurzem & 
Aarseth 1996). The reason is unclear, but see a discussion 
in subsections 3.2 and 3.3. In any system with even a slight 
external tidal field such binaries will be removed very fast. 

For runs with primordial binaries we discuss three dif- 
ferent cases. First, we use an artifical model with two- 
components, one consisting of single, the other of binary 



© 0000 RAS, MNRAS 000, 000-000 



20 M. Giersz and R. Spurzem 



stars, which only interact by relaxation with each other and 
with the single (all close 3b and 4b encounters were sup- 
pressed) . We check the rate of mass segregation and core col- 
lapse time by comparison with a continous two-component 
anisotropic gaseous model (Louis & Spurzem 1991, Spurzem 
1994), just in the same way as HA92 checked their A'^-body 
models by comparison with one of Reggie's gaseous models. 
We find that our stochastic Monte Carlo model of the bi- 
nary component provides a good match with the expectation 
from the gaseous model. While in Paper I we only checked 
the dynamical friction of one binary in a system of single 
stars, this is a more significant test of the self-consistent 
treatment of relaxation in our model. 

Second, a set of four models corresponding to the largest 
published full A''-body models with primordial binaries by 
Heggie & Aarseth (1992, HA92) was performed. We find 
a good agreement of our stochastic Monte Carlo model for 
most features, as they were published in HA92. The remain- 
ing differences regarding the binary number evolution, cen- 
tral potential, and fraction of hard binaries in the system 
between their and our models can be related to stochastic 
variations, as they usually occur in individual Monte Carlo 
realizations of the system. 

Third, and as a final goal of this paper, we have per- 
formed a self-consistent model of 30.000 binaries evolving 
internally and externally, surrounded by a live star cluster 
of 300.000 single stars, undergoing mass segregation, core 
collapse, and finally large amplitude gravothermal oscilla- 
tions. The binaries during all these evolutionary phases are 
subject to binary destructions and ejections by close 3b and 
4b encounters. This model is called Gao's run, to remind the 
reader of the pioneering work of Gao et al. (1991, GGCM91), 
which performed for the first time a self-consistent two- 
component Fokker-Planck model of binaries and single stars. 
Our aim was, in this paper, to redo this model as much as 
possible and reasonable, keeping most of the assumptions, 
such as the 3b and 4b interaction cross sections, the assump- 
tion of single mass (all binaries consist of mass components 
equal to the single stars), and the assumption of an isolated 
point-mass system. However, in addition to their model, we 
exploit in this framework the strengths of our stochastic 
Monte Carlo model, such that we are able to follow the in- 
dividual evolution of internal and external parameters of any 
binary in nearly the same detail as in an A'^-body model, pro- 
vide snapshots of the position and binding energy of all the 
binaries at crucial times of the evolution (first mass segre- 
gation, binary heating and destruction phase, gravothermal 
oscillation phase), which can be seen by the interested reader 

as a movie under 

ftp: //ftp. ari.uni-heidelberg.de/ 

pub/ spurzem/movie? .mpg 
where "?" refers to 1 or 2, depending whether one wants to 
see the movie in the style of Fig. ^ or ^ Also we do not use 
the rather artificial procedure of picking binaries in the en- 
ergy bins for close encounters as used by GGCM91, because 
in their model no better procedure could be found. Since we 
know positions of binaries at any time, we can select proper 
probabilities for close encounters of the 3b and 4b kind to 
occur and find a much slower destruction rate of soft and 
intermediate binaries than GGCM91. As a consequence the 
quasistationary binary burning phase extends over a time 
in our model (up to 150 frho ) which is about three times 



larger than the corresponding time in GGCM91. After that, 
gravothermal oscillations start, and at about 280 trho all pri- 
mordial binaries left the core and its vicinity. Only a fossil 
collection of some 400 binaries can be found in the very outer 
halo, as a final trace of the initial primordial binary popula- 
tion. During the maximum density phase new binaries start 
to be formed by 3b interactions. Of the initially 30.000 bi- 
naries, 15.000 were ejected by 3b and 4b encounters, and 
another 15.000 destroyed by 4b encounters. The single stars 
originating from such binary destructions are partly heat- 
ing the cluster, partly escaping, depending on the binding 
energy of the harder binary in such encounters. On average 
we find that for any binary, whether destroyed or ejected, 
about three single stars are ejected. Qualitatively we find 
the same phases as GGCM91, but our evolution time scales 
are a factor of 3 longer. Though the evolutionary time scale 
is much longer than the Hubble time, we think our models 
are still (already on the present idealized level) able to pro- 
vide some insight into the physical evolution of large star 
clusters with many primordial binaries, at least if accepting 
the point mass approximation. Note, that the inclusion of 
more realism (see below) according to previous experiences 
always shortens the evolutionary lifetimes of clusters, and 
stellar mass loss will change the cluster's and the individ- 
ual binaries evolution dramatically. Nevertheless we consider 
our work as an important and indispensable step towards 
such models with very large binary and single star numbers. 

Admittedly our model is still not realistic for a globular 
cluster in many respects. The evolution of hard binaries dur- 
ing close encounters is known to be very different from that 
of equal mass stars; stellar evolution and finite size efi'ects 
will dramatically alter the channels through which primor- 
dial binaries are processed (see Hut, McMillan & Romani 
1992). Tidal fields may strip the outer areas of the clus- 
ter, especially those where we find our remaining binaries at 
the end. Nevertheless we stress, that in our model we have 
achieved a breakthrough in the homogeneous, self-consistent 
treatment of very many binaries in a large star cluster. 
We model the cluster evolution during 0.8 million A'^-body 
time units, which are some 260.000 initial half-mass cross- 
ing times, without either using special hardware or other 
than the fundamental assumptions of spherical symmetry 
and dominance of small angle two-body encounters (which 
are very robust assumptions in large A'^ systems). Such a 
job is enormous, even for the next generation of Petafiop 
computers if done by a direct A'-body model. In this sense 
our model has proven its usefulness, uniqueness, and pio- 
neering capabilities, which gives credit to those people, who 
have been developing the Monte Carlo method (Henon 1971, 
Spitzer 1975, Stodolkiewicz 1982, 1986), and it shows the rel- 
evance of the new Monte Carlo models developed recently 
(Giersz 1996, 1998). 

The yet missing effects (multi-mass, stellar evolution, 
tidal fields, more detailed cross-sections, obtained by direct 
modelling) will be included into our models for the near 
future, and do not pose a fundamental challenge. But in or- 
der to do the proper astrophysical study, we rather prefer 
to do it step by step, connecting our models with Reggie's, 
Giersz's and Gao's runs, discuss their reliability, agreement 
and disagreement and the physical reasons of it, to under- 
stand what goes on in our model, as compared to just start 
with one big model containing everything. 
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